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Abstract 

We review recent experiments on dewetting thin films of evaporating colloidal nanoparticle suspensions 
(nanofluids) and discuss several theoretical approaches to describe the ongoing processes including coupled 
transport and phase changes. These approaches range from microscopic discrete stochastic theories to 
mesoscopic continuous deterministic descriptions. In particular, we focus on (i) a microscopic kinetic 
Monte Carlo model, (ii) a dynamical density functional theory and (iii) a hydrodynamic thin film model. 

Models (i) and (ii) are employed to discuss the formation of polygonal networks, spinodal and branched 
structures resulting from the dewetting of an ultrathin 'postcursor film' that remains behind a mesoscopic 
dewetting front. We highlight, in particular, the presence of a transverse instability in the evaporative 
dewetting front which results in highly branched fingering structures. The subtle interplay of decomposition 
in the film and contact line motion is discussed. 

Finally, we discuss a simple thin film model (iii) of the hydrodynamics on the mesoscale. We employ 
coupled evolution equations for the film thickness profile and mean particle concentration. The model is 
used to discuss the self-pinning and de-pinning of a contact line related to the 'coffee-stain' effect. 

In the course of the review we discuss the advantages and limitations of the different theories, as well as 
possible future developments and extensions. 

The paper is published in: /. Phys.-Cond. Mat. 21, 264016 (2009), 
in the Volume "Nanofluids on solid substrates" and can be obtained at 
http://dx.doi.org/10.1088/0953-8984/21/26/264016 
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I. INTRODUCTION 



The patterns formed in dewetting processes have attracted strong interest since Reiter analysed the 
process quantitatively in the early nineties. In these experiments, that proved to be a paradigm in 
our understanding of dewetting, a uniform thin film of polystyrene (tens of nanometers thick) is 
deposited on a flat silicon oxide substrate is brought above the glass transition temperature. The 
film ruptures in several places, forming holes which subsequently grow, competing for space. As a 
result, a random polygonal network of liquid rims emerges. The rims may further decay into lines 
of small drops due to a Rayleigh-type instability [HHH. The related problems of retracting contact 
lines on partially wetting substrates and the opening of single holes in rather thick films have also 
been studied [HH. 

Subsequent work has mainly focused on many different aspects of the dewetting process for simple 
non-volatile liquids and polymers (for reviews see Refs. [l6]-[8l). All stages of the dewetting of a 
film are studied: the initial film rupture via nucleation or a surface instability (called spinodal 
dewetting) [[ni9l-[T3l, the growth process of individual holes HIHUII, the evolution of the resulting 
hole pattern [[3l [131, and the stability of the individual dewetting fronts [fT7] - fT9l . We note in 
passing, that descriptions of dewetting patterns may also be found in historic papers, particularly 
for the dewetting of a liquid film on a liquid substrate. Tomlinson [[20l footnote 18 on p. 40] 
considered turpentine on water and Marangoni [|2T1 p. 352f] oil on water. 

More recently, interest has turned to the dewetting processes of solutions and suspensions. How- 
ever, these systems have not yet been investigated in any great depth. Such systems are compli- 
cated because their behaviour is determined by the interplay between the various solute (or colloid) 
and solvent transport processes. Furthermore, the solvents that are used often evaporate, i.e., one 
has to distinguish between 'normal' convective dewetting and evaporative dewetting. A number 
of experiments have been performed employing (colloidal) solutions of polymers [|22l - l25l . macro- 
molecules like collagen and DNA [|26] - l3ni and nanoparticles [|32] - l40l . The latter are sometimes 
referred to as 'nanofluids'. The initial focus of much of the research in the field has been on 
investigating the structures that are formed which are similar to the ones observed in the 'classi- 
cal' dewetting of non- volatile liquids. Labyrinthine structures and polygonal networks result from 
spinodal dewetting and heterogeneous nucleation and growth, respectively. They are 'decorated' 
with the solute and therefore conserve the transient dewetting pattern as a dried-in structure when 
all the solvent has evaporated [|28l l34l . The picture is, however, not complete. The solute may 
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also shift the spinodal and binodal lines as compared to the locations of these lines in the phase 
diagram for the pure solvent [41 J. As a consequence, the solute concentration influences the hole 
nucleation rate. More importantly, the solute particles may also destabilise the dewetting fronts. 
As a result, one may find strongly ramified structures in all three systems [|23l l25l l40l |42| . A 
selection of images exhibiting some of the possible structures is displayed in Fig|T] 
For volatile solvents, the contact lines retract even for wetting fluids. It has been found that such 
evaporatively receding contact lines may deposit very regular line or ring patterns parallel to the 
moving contact line [|24l l43ll . The deposition of a single ring of colloids from a evaporating 
drop of colloidal suspension is well known as the 'coffee stain effect' [|44l . Detailed investiga- 
tions reveal the emergence of rich structures including multiple irregular rings, networks, regular 
droplet patterns, sawtooth patterns, Sierpinski carpets, and - in the case of DNA - liquid crys- 
talline structures [l22l[30ll45l - |49l . The deposition of regularly spaced straight lines orthogonal to 
the moving contact line has also been reported [|50l . Droplet patterns may as well be created em- 
ploying solvent-induced dewetting of glassy polymer layers below the glass transition temperature 

mMi. 

Note that the dewetting of pure volatile liquids has also been studied experimentally [ |54| and 
theoretically Il55l458ll . In this case, different contact line instabilities have been observed for evap- 
orating liquid drops [|59ll60l . 

In the present article we review and preview the experiments and in particular the various mod- 
elling approaches for dewetting suspensions of (nano-)particles in volatile partially wetting sol- 
vents. After reviewing the basic experimental results in Section |Il} we discuss in Section III 



sev- 



eral theoretical approaches. In particular, we present a kinetic Monte Carlo model in Section III A 



a dynamic density functional theory in Section |IIIB[ and a thin film evolution equation in Sec- 
tion [lirC Finally, we conclude in Section |IV] by discussing advantages and shortcomings of the 
individual approaches and future challenges to all of them. 



n. EXPERIMENT WITH NANOPARTICLE SOLUTIONS 



We focus on experiments that use monodisperse colloidal suspensions of thiol-passivated gold 
nanoparticles in toluene [|33ll34ll37] - |40l[6T|| . The gold core of 2 - 3 nm diameter is coated by a layer 
of alkyl-thiol molecules. The length of the carbon backbone of the thiol used in the experiments 
ranges from 6 to 12 carbon atoms {Cq to C12) [l40ll . By varying the chain length, one can control 



4 




FIG. 1: (Colour online) Images of strongly ramified dewetting structures obtained using Atomic Force 
Microscopy in the case of (a) an aqueous collagen solution on graphite (courtesy of U. Thiele, M. Mertig 
and W. Pompe; see also Ref. f42l. Image size: 5/umx5/im); (b) poly(acrylic acid) in water spin-coated onto 
a polystyrene substrate (reprinted with permission of John Wiley & Sons, Inc. from Ref. ll23l : copyright 
John Wiley & Sons, Inc. 2002; Image size: 2.5/^mx2.5/xm); and in both (c) and (d), a solution of gold 
nanoparticles in toluene, spin-coated onto native oxide terminated silicon substrates (scale bars given in 
panels). In all the images the lighter areas correspond to the deposited solute and the dark areas to the 
empty substrate. 
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to a certain extent the particle-particle attraction. Normally, the solution is deposited on to a plain 
silicon substrate that is covered by the native oxide layer only ll34ll . However, one may locally 
change the wetting behaviour of the solvent by further oxidising the substrate (381 • By adding 
excess thiol one can also vary the properties of the solvent [|40ll . 

Two different procedures are employed for the deposition of the solution on to the substrate: spin- 
coating or a meniscus technique [|6T1 l62l . The choice is important as it strongly influences the 
evaporation rate and, as a result, the pattern formation process. When using spin-coating, one finds 
that directly after deposition, evaporation competes with dewetting until all the solvent has evapo- 
rated. The resulting deposits of nanoparticles are imaged by atomic force microscopy (AFM). For 
spin-coated films, the evaporation rate is high and structuring is normally finished before the spin- 
coater is stopped. Conversely, the solvent evaporation rate is strongly decreased when employing 
the meniscus technique [|6T1l . i.e., by depositing a drop of solution on a Teflon ring that is wetted by 
the solvent. This allows for a better control of the process and enables the use of contrast-enhanced 
microscopy to observe the dewetting process in situ POl . All pattern formation is confined to the 
region of the receding contact line of toluene, silicon and air. With both techniques one may find 
mono-modal or bi-modal polygonal networks [[34l, labyrinthine spinodal structures, or branched 
patterns (see Fig. [T]). The meniscus technique allows for the study of branched structures in a 
more controlled manner. The work in Ref. [40J indicates that fingering strongly depends on the 
interaction strength of the particles, i.e., on the chain length of the thiol molecules coating the gold 
cores. For short chains (C5 and Cg) no formation of branched structures is observed. At similar 
concentrations, well-developed branched structures are formed for longer chains (Cio and C12). 
For even longer chains (C14), however, one again finds less branching. It also depends on the 
amount of excess thiol in the solvent (for details see Ref. [40]). 

When following the evolution of the branched patterns in situ (see the complementary video 
material of Ref. ^4011 ). one clearly observes that different processes occur on different lenght 
scales. First, a macroscopic dewetting front recedes, leaving behind a seemingly dry substrate. 
The macroscopic front can be transversely unstable resulting in large-scale (> lOOyum) strongly 
anisotropic fingered structures. For fronts that move relatively quickly these macroscopic struc- 
tures cover all the available substrate. However, when at a later stage the macroscopic front be- 
comes slower, those fingers become scarce and 'macroscopic fingering' finally ceases. At this 
stage it is possible to appreciate that the seemingly dry region left behind by the front is not at all 
dry, but covered by an ultrathin 'postcursor' film that is itself unstable. The thickness of this film 
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is similar to the size of the nanoparticles. At a certain distance from the macroscopic front, the 
ultrathin film starts to evolve a locally isotropic pattern of holes. The holes themselves grow in an 
unstable manner resulting in an array of isotropically branched structures as shown, e.g., above in 
Fig. [T| This indicates that at least some of the patterns described in the literature may have arisen 
from processes in similar ultrathin 'postcursor' films. 

The existence of the ultrathin 'postcursor' film is an experimental finding that can be drawn on 
when choosing a theoretical approach to account for the pattern formation (see below). Note how- 
ever, that at the moment there exists no explanation for its existence. A possible hypothesis is 
that the substrate strongly attracts the nanoparticles. As a result they form a dense suspension 
layer having a thickness roughly equal to the diameter of the nanoparticles. The observed meso- 
scopic dewetting front then actually correspond to an autophobic dewetting of a low concentration 
suspension from the higher concentration suspension on the surface of the substrate. 



III. MODELLING APPROACHES 



Models of dewetting thin films of pure liquids or polymers are often based on thin film hydro- 
dynamics. Starting from the Stokes equations, together with continuity and boundary conditions 
at the substrate and free surface, one applies a long-wave approximation (assuming small surface 
slopes and contact angles) |l8l[63]| and obtains a non-linear evolution equation for the film thickness 
profile h{x, y, t). In the case of volatile liquids one finds Il55l - [58l[64| 
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with the mobility functions Qc{h) = h^/3r] > (assuming Poiseuille flow in the film and no slip 
at the substrate; t] is the dynamic viscosity) and Qc > for the convective and evaporative part 
of the dynamics, respectively. Qe is a rate constant that can be obtained from gas kinetic theory 
or from experiment [|57l . Note that Eq. ([T]) only applies if the pressure in the vapour above the 
film is close to the saturation pressure. For alternative expressions that are used to describe the 
non-conserved evaporative dynamics see, e.g., Refs. Il56l[57l I65l - l69l . Finally, V = {d^^dy), and 
dt, dx and dy denote partial derivatives w.r.t. time and the coordinates. 

Focusing on the influence of capillarity and wettability only, the energy functional F[h] is given 
by 



F[h] 



dx / dy 



1 



(Vhf + f{h) - fih 



(2) 



where 7 is the hquid-gas surface tension and f{h) is a local free energy term that describes the 
wettability of the surface. Since fi corresponds to a chemical potential, the term fih may either bias 
the system towards the liquid or towards the gas state. The variation of F w.r.t. h gives the pressure. 
It contains the curvature (Laplace) pressure —'jAh and the disjoining pressure Ii{h) = —dhf{h). 
Many different forms for the latter are in use (see, e.g., Refs. [|4l [81 16311701 - 1731 '). 
For the present system a thin film description using Eq. ([T]) is not appropriate because the nanopar- 
ticles are not taken into account. However, under certain conditions one can augment equation ([T]) 
for the evolution of the film thickness by coupling it to an equation for the evolution of the mean 
particle concentration. The resulting model is able to describe the behaviour of an evaporating so- 



lution on the meso- and macroscale. Such an approach is briefly discussed below in Section IIIC 



We should expect such a model to describe the mesoscopic dewetting front discussed above. How- 
ever, the theory is less suited to a description of the dewetting dynamics of the ultrathin postcursor 
film. 

The dewetting of the ultrathin film of highly concentrated suspension may be described by a dis- 
crete stochastic model such as, for instance, a kinetic Monte Carlo (KMC) model based solely on 
evaporation/condensation dynamics of the solvent and diffusion of the solute [l35l[39ll4Tl| . The va- 
lidity of this strong assumption regarding the relevant transport processes can be confirmed from 
an estimate based on Eq. ([T]): The pressure p = 5F/5h drives convection and evaporation. The 
convective mobility is proportional to h?, i.e., it is large for thick films but decreases strongly with 
reduced film thickness. The evaporative mobility, however, is a constant, implying that evapora- 
tion will dominate below a certain (cross-over) thickness. For the parameter values of Ref. Il57ll 
and a small contact angle (^ 0.01), the cross-over thickness is in the range of 1-5 nanometers. 
This estimate justifies the neglect of convective transport in a description of the postcursor film 
and may explain why one has such good agreement between the experimentally observed patterns 
and the patterns obtained from a purely two-dimensional (single layer) kinetic Monte Carlo model 
||35ll . We introduce the KMC model below in Section [mAl 

In several respects, however, the kinetic Monte Carlo model is rather simplistic, hmiting its po- 
tential applications. For instance, the thermodynamic chemical potential as well as any wetting 
interaction of the solvent with the substrate are collected in a single parameter - an effective chem- 
ical potential. This implies that any influence of a disjoining pressure is 'smeared out' over the 
whole system and that no distinction between the short- and the long-range parts of the disjoining 
pressure is possible. It is furthermore based on the assumption that evaporation/condensation is 
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the dominant dynamic process, but does not allow one to probe this assumption. In Section [Til B 



we show how one may develop a dynamical density functional theory (DDFT) that describes the 
system at a similar level to the KMC. However, the DDFT may also be easily extended to include 
other effects such as fluid diffusion, that the KMC does not incorporate. 



A. Kinetic Monte Carlo model 

The kinetic Monte Carlo model for two-dimensional dewetting nanofluids [|33l was first proposed 
in Ref. [|35l and extended to include next- nearest neighbour interactions in [[37l . The two key 
assumptions used are: (i) the relevant processes can be mapped on to a two-dimensional lattice 
gas model, thereby neglecting continuous changes in the thickness of the evaporating film, and (ii) 
all relevant dynamics results from diffusing nanoparticles and evaporating/condensing solvent. 
The model builds on an Ising-type model for the liquid-gas phase transition. The surface is divided 
up into a regular array of lattice sites whose size is dictated by the nanoparticles. One then con- 
siders each lattice site to be occupied either by a nanoparticle, liquid or vapour. This effectively 
maps the system onto a two-dimensional two-component lattice gas having two fields n and /. The 
resulting three possible states of a cell are: liquid (/ = l,n = 0), nanoparticle {I = 0,n = 1), 
and vapour (/ = 0, n = 0, i.e., cell empty). The energy of an overall configuration is given by the 
hamiltonian 

E = Yl ^^^3 - Y ^ ^'^^ ~ y ^ ^'^^ - f^J2^i 

<ij> <ij> <ij> i 

where J2<ij> denotes a sum over nearest neighbour pairs and su, Snn and are the liquid-liquid, 
particle-particle and liquid-particle interaction energies, respectively. Fixing the three interaction 
strength parameters Eu, £nn, ^ni and the effective chemical potential fi determines the equilibrium 
state of the system. We choose cu as unit of energy - i.e. we set Eu = 1. 

The hamiltonian determines the equilibrium state and the energy landscape of the system. How- 
ever, as the system 'dries in' during the course of the solvent evaporation, the final nanoparticle 
configurations do not necessarily represent equilibrium structures. This implies that the system 
dynamics is of paramount importance. It is determined by the possible Monte Carlo moves, their 
relative frequencies, and the probabilities for their acceptance. Two types of moves are allowed: (i) 
evaporation/condensation of liquid and (ii) diffusion of nanoparticles within the liquid. A mobility 
M corresponds to the ratio of cycles of particle and solvent moves and reflects the physical ratio of 
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time scales for evaporation and diffusion. A large mobility M indicates fast diffusion as compared 
to evaporation. A trial move is accepted with the probability pacc = min[l, exp{—AE/kT)] where 
k is the Boltzmann constant, T the temperature and AE is the change in energy resulting from the 
potential move. Note that particles are only allowed to move into wet areas of the substrate, i.e., 
onto cells with / = 1. This models zero diffusivity of the particles on a dry substrate. The replaced 
liquid fills the site left by the nanoparticle. 

Without nanoparticles, the behaviour of the model is well known as it reduces to the classical 
two-dimensional Ising model ^TM . For kT < kTc ~ 0.567 liquid and vapour coexist when fi = 
l^coex = —2. For yU > — 2 [yU < —2] eventually the liquid [vapour] dominates. A straight liquid- 
gas interface will recede [advance] for ji < —2[fi > —2], i.e. one finds evaporative dewetting 
[wetting] fronts. If one starts, however, with a substrate covered homogeneously by the liquid, 
for ji < —2 the film will dewet via a nucleation or spinodal-like process. If the nanoparticles are 
present, they form dried- in structures when all the liquid evaporates. The final structures do not 
normally change any further - at least on short time scales. However, if the liquid wets the particles 
(i.e. is attracted to the particles), over long times there might be a coarsening of the structures, 
facilitated by the adsorbed liquid. The dried-in patterns depend on the particular pathway taken by 
the evaporative dewetting process. They range from labyrinthine to polygonal network structures 
or holes in a dense particle layer. Some typical patterns are displayed in Fig. [2[ for cases when 
the average surface coverage of the nanoparticles p^" = 0.2. Panels (a) and (b) result from a 
spinodal-like and nucleation and growth process, respectively. At first sight they look very similar 
to the patterns seen for the pure solvent and one might argue that the particles solely act as passive 
tracers and preserve the transient volatile dewetting structures of the solvent. This was suggested 
in Refs. [l26] - l28ll for dewetting collagen solutions. However, panels (c) and (d) indicate that the 
particles may at times play a rather more significant role. When the diffusion of the particles is 
slow, the evaporative dewetting fronts become transversely unstable and may result in strongly 
ramified patterns. This instability is caused by the nanoparticles. The lower their mobility, the 
stronger the fingering effect, i.e., there are more fingers in (c) than in (d) because in the latter the 
mobility is larger. 

The front instability is intriguing as it results in strongly branched structures. As the dewetting 
front moves, new branches are continuously created and existing branches merge at the moving 
contact line. However, the mean finger number in the streamwise direction of the resulting ramified 
pattern is a constant. This behaviour is in contrast to the front instabilities found for dewetting 
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FIG. 2: Typical KMC results for the final dried-in nanoparticle structures resulting from the evaporative 
dewetting processes of nanoparticle solutions (nanofluids) in the case of (a) a spinodal-like process at ^ = 
—2.55, (b) nucleation and growth of holes at /i = —2.3, (c) unstable fronts at /i = —2.3 and low mobility 
M = 5, and (d) unstable fronts at = —2.3 and medium mobility M = 10. The starting configuration in 
(a) and (b) is a homogeneous liquid film with uniformly distributed particles whereas in (c) and (d) a hole 
at the center is nucleated 'by hand'. The remaining parameters are (a,b) M = 50, €ni = 2.0, e„„ = 1.5, 
p"^ = 0.2, kT = 0.3, MC steps= 500, domain size 1200 x 1200; (c,d) e„„ = 2.0, e„/ = 1.5, = 0.2, 
kT = 0.2, MC steps= 3000, domain size 1200 x 1200. Lattice sites occupied by particles are coloured 
black, and the empty sites are coloured white. 
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polymers which only result in fingers without side-branches [|75l or fields of droplets left behind 
JH. 

A quantitative analysis shows that the mean number of fingers depends only very weakly on the av- 
erage concentration of the nanoparticles p"^; only the mean finger width increases with increasing 
concentration. However, decreasing the mobility (i.e., decreasing the diffusivity of the particles) 
leads to a much denser finger pattern and also causes the front instability to appear at an earlier 
stage, i.e., when the front instability is in its initial linear regime, it has a higher growth rate and a 
smaller characteristic wavelength (cf. Fig.[2]^c) and (d)). Decreasing the effective chemical poten- 
tial (increasing its absolute value) has a similar but less strong effect. For details see [|4T]| . These 
findings lead to the conclusion that the determining factor for the front instability is the ratio of 
the time-scales of the different transport processes. In particular, the front becomes more unstable 
when the velocity of the dewetting front increases as compared to the mean diffusion velocity of 
the nanoparticles. 

If the particle diffusivity is low, the front 'collects' the particles, resulting in a build up of the 
particles at the front that itself is slowed down. This makes the front unstable and any fluctuation 
along the front will trigger a transverse instability that results in an evolving fingering pattern. This 
happens even when the particle-liquid and particle-particle attractive interactions do not favour 
clustering (i.e. demixing of the liquid and the nanoparticles). In this regime, the instability is a 
purely dynamic effect and energetics plays no role in determining the number of fingers. We call 
this the 'transport regime' . 

To illustrate the influence of energetics (characterized by the interaction parameters £ij) on finger- 
ing in Fig. [3] we display the dependence of the mean finger number on particle-liquid interaction 
strength Eni- For e„; > 1.5 the mean finger number < / > is nearly constant; this is the trans- 
port regime. However, on decreasing Sni below 1.5, we observe a marked increase in the value 
of < / >, indicating that energy plays an important role in determining the number of fingers in 
this regime. In this parameter range, demixing of particles and liquid occurs at the moving front 
and increases its transverse instability. In this 'demixing regime', the wavelength of the fingering 
instability is determined by the dynamics and the energetics of the system. Decreasing e„/ further 
(below 1.4 in Fig. [3]) one first observes in regime (iii) a slight decrease in the average finger num- 
ber. This is a geometric effect resulting from our one-dimensional finger counting routine: The 
fingers increasingly break up and the dried-in pattern looks progressively isotropic. In regime (iv), 
the measure (/) does not represent a finger number but instead indicates a decrease in the typical 
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distance between particle clusters resulting from the demixing process that occurs already in the 
bulk liquid and is not related to the front instability at all. Note that one finds a similar sequence 
of regimes (i) to (iv) when increasing the particle-particle interaction strengths for fixed Sni (see 
Ref. [41 J) for further details. 
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FIG. 3: (Colour online) Dependence of the mean finger number left behind by the unstable dewetting 
front on the particle-liquid interaction strength Sni- The regions marked (i) to (iv) are discussed in the 
main text. The insets display typical snapshots obtained in the four different regions. Particles are black, 
liquid is grey (green online) and the empty substrate is white. The remaining parameters are kT = 0.2, 
M = 20, ;u = -2.2, p^^ = 0.1, enn = 2.0, domain size 1200 x 1200. For the insets, from left to right, 
eni = 1.2,1.4,1.45,1.8. 



We note also that the fingering process may be viewed as self-optimising the front motion - i.e. 
the front keeps its average velocity constant by expelling particles into the fingers. A similar effect 
exists for dewetting polymer films [|T8l . where liquid is expelled from the growing moving rim 
which collects the dewetted polymer. There, the surplus liquid is left on the surface as a droplet 
pattern. 

The kinetic Monte Carlo model is a very useful tool that helps one to understand the pattern 
formation in drying nanoparticle suspensions. One has, however, to keep in mind the restrictions 
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on the model (see above). The purely two-dimensional character of the KMC was extended to 
a 'pseudo three-dimensional' one by making the effective chemical potential dependent on the 
mean liquid coverage [|38l . As the latter is related to a mean film thickness, this corresponds to 
the introduction of a 'global' thickness-dependent disjoining pressure into the evaporation term 
without an explicit consideration of a film thickness. The amended model can reproduce bimodal 
structures that are beyond the scope of the purely two-dimensional model [|38l [39ll . Fully three- 
dimensional models are also discussed in the literature [l76ll77l . 

B. Dynamical Density Functional theory 

The limitations of the kinetic Monte Carlo model introduced in the previous Section are related 
to its character as a two-dimensional lattice gas with only three states: gas, liquid or particle. 
This implies that (i) no liquid can be transported to a site on the surface already filled with liquid, 
i.e., diffusion of the liquid can not be incorporated in a sensible way and (ii) one is not able to 
distinguish between the influence of the short- and the long-range parts of the interactions with the 
substrate, as all such interactions are absorbed into the effective chemical potential. 
However, using dynamical density functional theory (DDFT) [iTSl - ISSl one can develop a model 
for the processes in the ultrathin postcursor film without these limitations, although here we limit 
ourselves to developing the theory at the level of the KMC and solely discuss how to extend it to 
incorporate the influence of the liquid diffusion over the surface. Such a DDFT model describes 
the coupled dynamics of the density fields of the liquid pi and the nanoparticles p„. The densities 
pi and pn are defined as the probabilities of finding a given lattice site on the surface to be occupied 
by a film of liquid or by a nanoparticle, respectively. Note that the probability densities correspond 
to number densities as we use the lattice spacing a = 1 as our unit of length. 
To develop the DDFT, one must first derive the underlying free energy functional F[p/,p„], and 
secondly, devise dynamical equations for both density fields that account for the conserved and the 
non-conserved aspects of their dynamics, i.e., transport and phase change processes, respectively. 
For a system governed by the hamiltonian (|3]), we may construct a mean-field (Bragg-Williams) 
approximation for the free energy of the system [1781 [841 which contains an entropic contribution 
and contributions from the interactions between the different species (nanoparticles and liquid). 
The free energy is a semi-grand free energy, since the liquid is treated grand canonically (it is 
coupled to a reservoir with chemical potential p), whereas the nanoparticles are treated in the 
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canonical ensemble. The free energy functional is first defined on the original KMC lattice. How- 
ever, after re- writing the interaction terms employing gradient operators [TTSl one finally obtains 
the free energy functional for a continuous system 



F[puPn] = y"dr[/(p,,p„) + ^(VpO' + ^(Vp„y + £„/(Vp„)-(VpO-m 



(4) 



where 



fipuPn) = A;r[p,lnp/ + (1 - p;)ln(l - Pi)] 
+ kT[pn lnp„ + (1 - pn) ln(l - p„)] 

- 2eiipf - 2ennPn - 4:£nlPnPl- (5) 

Since the liquid may evaporate from the surface into the vapour above the surface, p is the (true) 
chemical potential of this reservoir and determines the rate of evaporation [condensation] from 
[to] the surface. Note that normally a free energy of the form in Eq. (|4]) is obtained by making a 
gradient expansion of the free energy functional of a continuous system ll84l . However, here we 
have made the mapping from the free energy of the lattice KMC system. 

The chemical potential for the nanoparticles may be determined from the functional derivative 
Pn = SF[pn, pi]/6pn{r). lu equilibrium it is constant throughout the system, but it may vary 
spatially in a non-equilibrium system, i.e., p„ = p„(r,t). We assume that the dynamics of the 
nanoparticles is governed by the thermodynamic force Vp„ - i.e. that the nanoparticle current 
is j = — M„p„Vp„5 where M„(p;) is a mobility coefficient that depends on the local density of 
the liquid. Combining this expression for the current with the continuity equation, we obtain the 
following evolution equation for the nanoparticle density profile 



dPn^ 
dt 



(6) 



SPn 

Note that this equation of motion may also be obtained by assuming that the nanoparticles have 
over-damped stochastic equations of motion [ [80] - i83l . Here, we assume that Mn{pi) = aQs{pi — 
0.5), where Qs{x) is a continuous function that switches smoothly from the value to the value 
1 at X = (i.e. it is essentially a smooth analogue of the Heaviside function). This ensures that 
the nanoparticles are immobile when the local liquid density is small (dry substrate) and have a 
mobility coefficient a when pi is high (wet substrate). 

For the evolution of the liquid density distribution we assume that the liquid is able to evaporate 
from the surface into the vapour (reservoir) above the surface (non-conserved dynamics) and may 
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FIG. 4: (Colour online) Density profiles for the situation where the substrate is covered by nanoparticles 
with average density p'^ = 0.3. The top row are the nanoparticle density profiles and the bottom row are 
the corresponding liquid density profiles at the times t/t; = 8 (left) and 80 (right), where = 
The parameters are kT/su = 0.8, Eni/eu = 0.6, e„„ = 0, a = 0.4MfV^, Mf = 0, pi{t = 0) = 0.9 ± ^ 
(where ^ represents white noise of amplitude 0.05) and {fi — fj.coex)/kT = —0.88, where the liquid exhibits 
spinodal decomposition-evaporation. 

also diffuse over the substrate (conserved dynamics). The conserved part is treated along the lines 
developed above for the nanoparticles. For the non-conserved part we assume a standard form 
[|85ll . i.e., the change in time of pi is proportional to — (/isurf(i", t) — /i) = —6F[pn,pi]/6pi{r) 
where /Xsurf (r, t) is the local chemical potential of the liquid at the point r on the surface at time t. 
This gives the evolution equation for the liquid density 



dpi 
dt 



V- 



M[piV 



SF[pn,pl] 

5pi 



SF[pn,pi] 



hi 



(7) 



where we assume that the coefficients Mf and are constants. 
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FIG. 5: (Colour online) Density profiles for the situation where the substrate is covered by nanoparticles 
with average density p"^ = 0.3 and with the liquid excluded from the region y < 0. The top row shows 
the nanoparticle density profiles and bottom row the corresponding liquid density profiles at the times 
t/ti = 1000 (left), 10000 (middle) and 30000 (right), where ti = l/kTMf'^a'^. The parameters are 
kT/eu = 0.8, Eni/eu = 0.6, e„„ = 0, a = 0.2MfV^, Mf = 0, pi{t = 0) = 0.9 ± ^ (where ^ represents 
white noise of amplitude 0.05) and {p — pcoexj/kT = —0.78. 

This theory allows us to study the time evolution of the evaporating film of nanoparticle suspension 
without some of the restrictions of the kinetic Monte Carlo model. Here, however, we illustrate its 
application in similar parameter regimes as used above for the KMC. We focus on two examples: 
(i) the spinodal dewetting of a initially flat film of nanoparticle suspension characterised by con- 
stant pi and pn (Fig.|4]); and (ii) the retraction of a dewetting front that is unstable with respect to 
a fingering instability (Fig.|5]). 

Fig.|4]presents two pairs of snapshots from a purely evaporative dewetting process deep inside the 
parameter region of the phase diagram where spinodal dewetting occurs. For small times the film 
becomes unstable showing a typical spinodal labyrinthine pattern with a typical wavelength. The 
nanoparticles concentrate where the remaining liquid is situated. However, they are 'slow' in their 
reaction: when pi already takes values in the range 0.08 - 0.83, the nanoparticle concentration 
has only deviated by about 25% from its initial value. The film thins strongly forming many 
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small holes. The competition for space results in a fine-meshed polygonal network of nanoparticle 
deposits. The concentration of particles is much higher at the network nodes - an effect that can 
not been seen within the KMC model. As the particles attract the liquid there remains some liquid 
on the substrate where the nanoparticles are. 

Fig. |5] gives snapshots of the evolution of a fingering instability for a retracting dewetting front. 
At early times the straight front shows a rather short-wave instability, about 16 wiggles can be 
seen. However, they are only a transient: the finger pattern coarsens rapidly till only about 7 
fingers remain. The fingering then becomes stationary, i.e., just as in the KMC, the mean finger 
number remains constant, although new branches are continuously created and old branches join 
each other. In general, the results on fingering agree well with results obtained using the KMC 
model pT|. From this we conclude that jamming of discrete particles is not a necessary factor 
for causing the instability, since the fingering is seen here in a continuum model with a diffusion 
constant that is independent of the nanoparticle concentration. The DDFT is better suited than the 
KMC for investigations of the early instability stages: they are more easy to discern without the 
discrete background noise of the KMC. Furthermore, one may perform a linear stability analysis of 
the one-dimensional undisturbed streamwise front profiles with respect to transverse perturbations 
(in analogy to the approach used in Refs. [fT9l[86l[87l ). 

C. Thin film hydrodynamics 

The previous two sections focused on two approaches to describe the experimentally observed 
patterning dynamics in the ultrathin postcursor film left behind by a mesoscopic receding dewet- 
ting front. Although both the kinetic Monte Carlo model and the dynamical density functional 
theory are able to describe well the processes in the ultrathin film, they can not be employed to 
describe mesoscale hydrodynamics. A relatively simple model for the latter can be derived in the 
framework of a long- wave or lubrication equation [[8l[63]|. We will illustrate here the approach 
by considering an isothermal situation where the nanoparticles are not surface active, i.e., they do 
not act as surfactants. For a model incorporating the effects of latent heat generation and surface- 
active particles resulting in thermal and solutal Marangoni stresses, see Ref. [|88ll . A description of 
spreading particle solutions incorporating a structural disjoining pressure has also been considered 
ll89ll . For related work on particle-laden film flow on an incline see Refs. ll90ll9T1l . 
One starts from the Stokes equations, together with continuity, no-slip boundary conditions at the 
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substrate and force equilibria at the free surface, and applies a long-wave approximation. Under 
the assumption that concentrations equilibrate rapidly over the film thickness, we obtain coupled 
non-linear evolution equations for the film thickness profile h(x,t) and the amount of nanoparticles 
per unit length hp = (ph, where is the volume concentration of the nanoparticles. Note, that hp 
corresponds to the local thickness of the nanoparticle layer when all the solvent is evaporated. The 
resulting evolution equation for the film thickness is Eq. ([T]) above and focusing on the influence 
of particle-independent capillarity and wettability only, the energy functional F[h] is given by 
Eq. ([2]) above. Note that the viscosity r] depends on the particle concentration. Following Refs. 
[|88l[89l|9D|92l we use the Quemada law for dense suspensions [|93ti95]| 

77(0) =Vo(l- ^) (8) 



where (pc = 0.64 corresponds to random close packing of spherical particles. For the nanoparticle 
volume per length hp = cph one obtains the following evolution equation: 

.5F' 



dMh) = V ■ 



cV- 

' 5h 



+ V ■ [D{cp)hV<p\ , (9) 



where the particle concentration dependent diffusion coefficient D{(j)) is related to the viscosity by 
the Einstein relation -D(0) = kT/67rRr]((f)), where R is the radius of the nanoparticles [|96ll . 
We illustrate results obtained employing this thin film theory using the single example of a re- 
ceding dewetting front for a partially wetting film. We use the disjoining pressure and material 
constants for the liquid considered in Ref. [|57l . where the evaporative and convective dewetting 
of a film of volatile liquid is studied. We add, however, the nanoparticles to the system. The 
expression that we employ for the local free energy term in Eq. ^ is: 

f{h)= ^2 +^pexpl 1, (10) 

where the parameters characterising the interaction between the liquid film and the surface are 
the apolar and polar spreading coefficients Slw and Sp, respectively, the Debye length Iq and the 
Bom repulsion length do [57J. The resulting disjoining pressure 11 = —dhf{h) allows for a stable 
precursor film (thickness /iprecursor) and also has a second (larger) thickness (ho) that corresponds 
to a secondary minimum of the underlying energy functional. See Refs. [[TTl l97l for studies of 
film and drop states for similar disjoining pressures. Our results are calculated for a system where 
the profiles only vary in one Cartesian direction (x), corresponding to a straight dewetting front. 
However, our results may also be interpreted as applying to a circular flat drop whose front remains 



19 




FIG. 6: Profiles of the final dried-in nanoparticle layer for the dewetting of a suspension of nanoparticles 
in a volatile solvent that partially wets the substrate for (a) high (Tl = 10~^), (b) medium (Q = 2 x 10~^) 
and (c) low (0 = 0.78 x 10~®) evaporation rates, for the case when x = H/Iq = 1.09, the lateral length 
scale is ^ = \/^/nH with k = {Sp/lo) exp{do/lo)H being an energy scale related to wettabiUty and the 
vertical length scale is H = y/2SLw / f^d^. The remaining dimensionless parameters are the evaporation 
number O = Qer]o£^/H^, the diffusion number T = D{0)r]o/HK = 10~^ and the dimensionless chemical 
potential M = H^/k = —0.0035. The system size is L = 19500£. Film thickness and hp in the plots are 
scaled by the precursor film thickness. 

circular throughout the dewetting and evaporation process. In this case one should interprete the 
coordinate x as the distance from the centre of the circular film. 

We start with a film of height ho of finite length sitting on a precursor film and assume that the film 
contains nanoparticles at constant concentration 0o- The chosen parameter values ensure that the 
film of thickness ho is linearly stable. As we do not incorporate noise, no nucleation of additional 
holes can occur (even with noise the probability would be extremely low). Without evaporation the 
film dewets 'classically' by a retraction of the initially step-like front. After a short time, surface 
tension smoothes the profile of the receding front and a capillary rim forms that collects all the 
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dewetted liquid. The front recedes until all liquid is collected in a central drop. Since no liquid 
evaporates [Q^c = in Eq. ([1])], the particle concentration does not change during the process. 
The situation changes when allowing for evaporation (Qnc > 0). Now the front may retract 
by convection and/or evaporation. Evaporation leads to the possibility of a strong increase in 
the particle concentration at the contact line as evaporation is strongest there. Due to the strong 
nonlinear dependence of the viscosity on the particle concentration, this may lead to a dramatic 
decrease of the convective contribution to the front velocity. For moderate evaporation rates, this 
may result in a (temporary) self-pinning of the front. Within the present basic model, the process 
can (after complete dry-in) result in three different basic deposition patterns: (i) for very fast 
evaporation rates, all other processes occur over time scales that are much larger. In particular, the 
effects of convective redistribution of the liquid are neglectable. As a result one finds that a nearly 
homogeneous film of nanoparticles of thickness hp = (j)oho is deposited (see Fig. [6]; a)). Convection 
only results in the small heap of material visible at the left hand side of Fig. [6]^a). The decrease 
in hp on the right side of Fig. [6]^a) arises due to the diffusion of particles to the right of the initial 
front position; (ii) for very low evaporation rates, the film dynamics is dominated by convective 
dewetting as this process acts on a much shorter time scale than evaporation. As a result, all the 
liquid is collected into a drop before evaporation slowly removes the remaining solvent. Under 
these conditions most of the nanoparticles are deposited in a single heap (see Fig.[6]^c)). Depending 
on the diffusivity, the heap might be highest at the centre or show a depression there; (iii) at 
intermediate evaporation rates, one may observe the deposition of a nanoparticle ring around a 
region with a nanoparticle film of much lower height. At the centre deposition might increase 
again (see Fig.[6]^b)). 

The most intriguing feature is the ring formation that has been observed experimentally for sus- 
pensions of very different particle sizes ranging from nanometers [|32l [36l l46l l47ll to hundreds of 
micrometers. Pinning of the contact line and thermal Marangoni effects are often mentioned as 
necessary conditions for the ring formation. The contact line pinning is often assumed to result 
from substrate heterogeneities. Film height and concentration profiles at various instants during 
the dewetting process are displayed in Fig.|7j The profiles are from before, at and after self-pinning 
of the contact line. In Fig. [8] we display a space-time plot for the complete process. At first, the 
front recedes in the same manner as when there is no evaporation, but now driven by convection 
and evaporation. A small capillary rim forms that collects all the dewetted liquid that does not 
evaporate. The particle concentration slowly increases at the contact line (Fig. ^a) and regime 
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FIG. 7: (Colour online) A sequence of profiles during a dewetting process with competing evaporation and 
convection that leads to the dried-in ring structure of nanoparticles displayed in Fig.|6jb). Profiles are at (a) 
before pinning (t = O.OST), (b) at self-pinning (t = O.IST), and (c) after depinning (t = 0.29T), where 
T = 3 X W^^T with r = tjo^H/k^ (T is of order of Is). The film thickness profiles h are the bold solid 
lines, the nanopaiticle concentrations cj) are the dotted lines and the nanoparticle layer height hp = h(j) are 
the dashed lines. The remaining parameters and scalings are as in Fig.[6];b). 

(i) in Fig. [8]). The concentration increases further and when it approaches random close packing 
0c, the viscosity diverges and the front pins itself. When pinned, further retraction only occurs 
through evaporation (Fig. [TJ^b) and regime (ii) in Fig. [8]). The front eventually depins and starts 
to move again, leaving a nanoparticle ring behind (Fig.|7]^c) and regime (iii) in Fig. [8]). However, 
the velocity is not as large as at the beginning, owing to the fact that the mean concentration of 
particles has increased. The remaining particles are transported to the centre and are deposited 
there when the remaining solvent evaporates (regime (iv) in Fig. [8]). 

The simple model used here shows, (i) that the contact line stops due to self-pinning by the de- 
posited particles and (ii) the Marangoni effect is not necessary for the ring formation. The model 
can easily be refined to account for solutal and/or thermal Marangoni effects [88J but self-pinning 
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FIG. 8: (Colour online) Space-time plots are given for (left) the film thickness h and (right) the nanoparticle 
layer height hp = hcf). The plot corresponds to the complete evolution resulting in the ring profile of 
Fig. [6jb). In both panels bright [dark] parts denote high [low] regions. The prominent central dark-bright 
border in the left panel indicates the change of the position of the contact line in time. Over time, four 
regimes can be distinguished: (i) fast motion before pinning, (ii) nearly no front motion during self-pinning, 
(iii) slow motion after depinning, and (iv) final evaporation from the center. 

should also be investigated further in the simple case presented here. 
IV. CONCLUSION 

We have discussed recent work on pattern formation processes in films and drops of evaporating 
suspensions/solutions of polymers and particles. After reviewing experiments on suspensions of 
thiol-coated gold nanoparticles in toluene we have focused on the modelling of the transport and 
phase change processes involved. A theoretical approach to the modelling of the hydrodynamics 
on the mesoscale has been described as well as more microscopic models for the dynamics in the 
observed nanoscopic 'postcursor' film. In particular, we have introduced (i) a microscopic kinetic 
Monte Carlo model, (ii) a dynamical density functional theory and (iii) a hydrodynamic thin film 
model. 

The kinetic Monte Carlo model and the dynamical density functional theory can both be used to 
investigate and understand the formation of polygonal networks, spinodal and branched structures 
resulting from the dewetting of an ultrathin 'postcursor' film that remains behind the mesoscopic 
dewetting front. They are, however, not capable of describing the dynamical processes in a meso- 
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scopic film. We have seen that the KMC model is able to describe the interplay of solute diffusion 
within the solvent and solvent evaporation/condensation. It also takes the liquid-liquid, liquid- 
particle and particle -particle interactions into account and therefore allows us to distinguish differ- 
ent regimes of the transverse (fingering) instability of the evaporative dewetting front: a transport 
regime where the instability is almost completely independent of the interaction strengths and 
a demixing regime where particles and liquid demix at the receding front thereby increasing its 
transverse instability. 

The dynamical density functional theory describes the coupled dynamics of the density fields of 
the liquid and the nanoparticles. In the form described above (i.e. based on the two-dimensional 
hamiltonian ([3])) we obtain a simple theory that allows us to study the time evolution of the evapo- 
rating ultrathin film and also to investigate the influence of processes such as surface diffusion by 
the liquid, which are not incorporated in the KMC model. However, it is straightforward to extend 
the theory to consider a fully three-dimensional fluid film, in which one can distinguish between 
short- and long-range interactions of solvent and/or solute with the substrate. We have, however, 
restricted the examples given here to situations that can also be described using the KMC model. 
A further exploration will be presented elsewhere. 

Finally, we have discussed a simple thin film model for the hydrodynamics on the mesoscale. It 
results from a long-wave approximation and consists of coupled evolution equations for the film 
thickness profile and the mean particle concentration. It has been used to discuss the self-pinning 
of receding contact lines that is related to the formation of rings of dried-in particles (coffee- 
stain effect) that frequently occurs when films or drops of solutions or suspensions dewet by the 
combined effects of convection and evaporation. 

One of the primary goals of researchers in this field, is the search for simple-to-use techniques 
that allow one to produce hierarchically structured functional layers for a wide range of applica- 
tions such as, e.g., organic solar cells [|98l . This means that the experiments advance very rapidly 
towards increasingly complex systems. For example, there have been investigations of the influ- 
ence of the phase behaviour on the drying of droplets of a suspension of hard-sphere colloidal 
particles and non-adsorbing polymer [|99l . of the instabilities and the formation of drops in evap- 
orating thin films of binary solutions HI 0011 that may lead to treelike patterns HlOlll . of effects of 
a secondary phase separation on evaporation-induced pattern formation in polymer films [I102L 
and of the influence of an imposed flow on decomposition and deposition processes in a sliding 
ridge of evaporating solution of a binary polymer mixture [I103II and of the influence of rather 
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fast evaporation [|104[ I105II . These complex experimental systems all represent systems of high 
practical interest that the theories presented here are not (yet) able to describe. Such experiments 
do, however, provide a strong motivation for further work to extend the theories presented here, as 
well as to develop new approaches. 

Let us finally mention that several topics were entirely excluded from our discussion here. First, we 
focused on a limited range of descriptions and did, for instance, not mention lattice Boltzmann, 
molecular dynamics or dissipative particle dynamics approaches that may also be employed to 
describe fluid suspensions [I106I - I1091 . Second, we have only discussed spatially homogeneous 
substrates. Patterned substrates are widely used in dewetting experiments [|38llllOl4ll2ll . Theoret- 
ical descriptions are well developed for the dewetting of films of pure non- volatile liquids on such 
substrates [l68l I1131 - I119L However, in the case of volatile liquids on heterogeneous substrates, 
much less work has been done. A third topic that we did not touch upon are possible continuum 
thin film approaches to demixing dewetting suspensions. We believe it is feasible to extend the 
diffuse interface theories such as model-H HI 2011 to include the influence of evaporation in dewet- 
ting nanoparticle suspensions. For instance, such models have already been adapted to describe 
demixing free surface films of polymer blends [|12H - I123]| . 
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